-
Notifications
You must be signed in to change notification settings - Fork 5
Support wraparound indexing in longitude #47
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: main
Are you sure you want to change the base?
Conversation
b57caba
to
fb36168
Compare
fb36168
to
d42f267
Compare
I think this works with the tests that you added here because the affine transformation is linear (symmetrical) and also because how numpy works with negative indices. I don't think this is will always work without something to handle wraparound, though. |
indexed = ds.pipe(assign_index) | ||
# We lose existing attrs when calling ``assign_index``. | ||
assert_equal(ds.sel(lon=[220, 240]), indexed.sel(lon=[-140, -120])) | ||
assert_equal(ds.sel(lon=220), indexed.sel(lon=-140)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
For example
assert_equal(ds.sel(lon=[0, 0]), indexed.sel(lon=[0, 359.9]))
should pass with some proper wraparound but raises
IndexError: index 180 is out of bounds for axis 2 with size 180
for indexed.sel(lon=359.9)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
hmmm.. using method="nearest"
by default is quite wrong, but CoordinateTransformIndex requires it
I would be nice to also have an option for shifting the range [0, 360] from/to [-180, 180]. |
d42f267
to
f4cd46a
Compare
Initially I thought that is what this PR was addressing. I guess I didn't realize that the affine "automatically" handles out-of-range values (for certain cases if I understand @benbovy's comment above correctly)! I was expecting a function or keyword argument to
I guess CRS is needed to check that coordinates are lon/lat and raise if the CRS is projected (where there is no wrapping)? https://pyproj4.github.io/pyproj/stable/api/crs/crs.html#pyproj.crs.CRS.is_projected In my mind there are two scenarios: 1. Selection using an affine that generates coordinates in [0,360] & user supplies out-of-bounds longitude that needs to be wrapped into the [0,360] range def wrap_coord(lon):
# same as converting from [-180,180] to [0,360] range
return lon % 360 da.sel(lon=-120) == ds.isel(lon=240)
da.sel(lon=365) == ds.isel(lon=5) 2. Selection on affine that generates coordinates in [-180,180] & user supplies longitude in [0, 360] range that needs to be converted da.sel(lon=240) == da.sel(lon=-120) def convert_to_negative_west(lon):
return ((lon - 180) % 360) - 180 |
Naive question: is it common to have unprojected lat/lon raster datasets with a rectilinear / no-rotation affine transform? It is pretty common for model outputs but they don't use affine transforms. |
FWIW Apache SIS has a special WraparoundTransform class. I guess one option could be to provide for convenience a similar class that would interoperate with |
Yeah, this is a lot more subtle than I initially realized, I won't be able to do it before SciPy. I am going to issue an alpha release now. |
Closes #26
Apparently we inherit this from
affine
, so no CRS needed. Does anyone understand this?